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The preparation of a quantum degenerate gas of heteronuclear molecules has been an outstanding 
challenge. We use path integral Quantum Monte Carlo simulations to understand the role of inter¬ 
actions and finite temperature effects in the protocol currently employed to adiabatically prepare a 
low-entropy gas of polar molecules in a lattice starting from an ultracold Bose-Fermi mixture. We 
find that interspecies interactions affect the final temperature of the mixture after the adiabatic 
loading procedure and detrimentally limit the molecular peak filling. Our conclusions are in agree¬ 
ment with recent experimental measurements [1] and therefore are of immediate relevance for the 
myriad experiments that aim to form molecules from dual-species atomic gases. 


Introduction. Polar molecules, interacting via long- 
range and anisotropic dipolar interactions, hold great 
promise as quantum simulators hosting exotic quantum 
phases [2, 3], as well as a diverse range of phenomena, 
ranging from quantum magnetism [4], to many-body lo¬ 
calization [5], to synthetic spin-orbit coupling [6, 7]. A 
necessary ingredient for simulating these these behaviors 
is to the ability to reach low entropy conditions. How¬ 
ever, despite the rapid experimental progress, a reliable 
method to form such a state has remained out of reach [8- 
10]. The two main obstacles to realizing a low-entropy 
state are the inapplicability of standard atomic cooling 
techniques to polar molecules and the requirement to 
suppress chemical reactions. A proposed solution to these 
problems is to form the molecules directly in a deep op¬ 
tical lattice through association after optimally loading 
a degenerate atomic gas mixture [11-20]. This protocol 
is shown schematically in Fig. 1. However, in order for 
this scheme to be a reliable pathway, it is necessary to 
have a high atom/molecule conversion efficiency by creat¬ 
ing a large region where the densities of the two atomic 
species overlap and correspond to exactly one atom of 
each species per lattice site. Temperature, interactions, 
and loading conditions can significantly limit the achieve¬ 
ment of this requirement. 

In this Letter, we use path integral Monte Carlo sim¬ 
ulations (QMC) based on the worm algorithm [21], as 
well as its two-worm extension, to investigate the effects 
of finite temperature and interspecies interactions dur¬ 
ing adiabatic loading from a dipole trap into an optical 
lattice. We show that the final temperature of the lat¬ 
tice system following adiabatic loading depends strongly 
on the strength and sign of interspecies interactions, and 
can be far from the ideal, zero-temperature regime. In 
contrast to previous studies, which do not treat the adi¬ 
abatic loading procedure and find that attractive inter¬ 
species interactions enhance the on-site densities [22, 23], 
our analysis predicts that both attractive and repulsive 
interactions can lead to substantial depletion of on-site 
densities of the bosonic species and takes into account the 
contributions of dimensionality, effective mass imbalance, 



FIG. 1. (Color online) Production of ground state molecules 
from a dual-species mixture. (Top line) A finite-temperature 
dual-species gas of bosonic species A atoms (red) and 
fermionic species B atoms (blue), is adiabatically loaded from 
a harmonic trap into an optical lattice in the presence of inter¬ 
species interactions. The initial temperature and interactions 
determine the temperature, density overlap, and peak filling 
of the mixture. (Bottom line) Following STIRAP, only those 
sites with exactly one atom of each species are converted to 
a ground state molecule, with some probability to excite the 
molecule to a higher lattice band (schematic wavefunctions in 
green box). 

and quantum statistics [24-26] . Additionally, we discuss 
the impact of imperfections in stimulated Raman Adia¬ 
batic Passage (STIRAP), which is used to convert Fesh- 
bach molecules to ground state molecules, on the ground 
state molecule production [11, 27]. We study the adia- 
baticity of this protocol, and provide approximate ana¬ 
lytical formulas which can be used to determine the prob¬ 
ability of promotion of molecules to higher bands during 
the STIRAP procedure. Since molecules in higher bands 
have a much larger tunneling rate, an appreciable higher 
band population can greatly impact dynamics involving 
molecular motion. We find that for experimentally rele¬ 
vant parameters there is no significant population trans¬ 
fer, provided the Lamb-Dicke criterion is satisfied. Our 
conclusions are corroborated by the recent experimen¬ 
tal observations reported in [1] and thus of fundamental 
importance for ongoing experimental efforts to achieve a 
high-filling lattice gas of ground state polar molecules. 

Single-component gases in the lattice We begin by 
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studying trapped single-species gases of bosons or 
fermions. We assume the bosons are trapped in a deep 
lattice where they can be described by the Bose-Hubbard 
Hamiltonian with an additional external harmonic con¬ 
finement: 

H = - ^ JA(aJ A aj,A + h.c.) + — ^ n itA { n i,A ~ 1) 
(ij) i 

-^/U,A™i,A- (1) 


Here, aJ A (a^ a) is the bosonic creation (annihilation) 

operator for species “A” and a = a \ A a i, a- The 
first and second terms in the Hamiltonian Eq. (1) are 
the tunneling J A and the on-site interaction, assumed 
repulsive with strength U >0, respectively. Finally, 
/U,a = M — Y,£=x,y,z w ^h where w% and & are the 
strength of harmonic confinement and the coordinate of 
site i along axis £, respectively, and fi is the chemical 
potential. 

In order to facilitate comparisons with the experi¬ 
mental data in [1], we will use 87 Rb for species A, but 
we stress that our conclusions apply for generic bosonic 
species described by the Bose-Hubbard model. The ef¬ 
fective external harmonic confinement includes the ad¬ 
ditional confinement due to the curvature of the lattice 
beams. The A atoms are assumed to experience a lattice 
potential Vo = 20£U, where Ea is the recoil energy of 
the A atoms. 

To characterize the density as a function of atom num¬ 
ber, we first match the entropy of the gas in the dipole 
trap to the one in the lattice and find the final temper¬ 
ature Tfi n /JA. The bottom right inset of Fig. 2 shows 
an example of this procedure for 3000 atoms. The solid 
green line shows the entropy of a weakly interacting Bose 
gas, which describes the A atoms in the harmonic trap. 
Blue triangles denote the entropy of the A atoms after 
loading into the lattice, given by 

(2) 

where E(T) is the system energy at temperature T. Both 
E(T) and E(0) are directly measured in our QMC simu¬ 
lations. An example of the entropy matching procedure 
for an initial temperature = 0.4T C , T c the critical tem¬ 
perature for Bose-Einstein condensation, is shown with 
dashed arrows. 

Next we determine the peak filling, f p . To compare 
with the experimental procedure of extracting the peak 
filling, we first integrate along z and mimic the experi¬ 
mental imaging resolution by applying a Gaussian filter 
with a width of 4 sites to the integrated density. We 
fit the resulting density to a Thomas-Fermi (TF) distri- 

bution n(x,y) = § f p a z [l - ( x/a x ) 2 - (y/cr y ) 2 ] 3/2 . The 
extracted f v values are plotted in Fig. 2 as a function of 



FIG. 2. (Color online) Peak filling of species A vs. atom 
number. (Main panel) The peak filling f p as a function of 
number of 87 Rb atoms, used as an example of bosonic species 
A, with a lattice depth of Vo = 20 Fa given for harmonic 
confinements of cj r (cj z ) — 27r40(260) Hz (filled triangles) 
and 27T35(227) Hz (filled circles). Empty diamonds are ex¬ 
perimental data [1]. The dashed (solid) line shows the on¬ 
site density at the center of the trap at T/J A — 15 for the 
weaker (stronger) harmonic confinement. (Lower inset) En¬ 
tropy matching procedure for 87 Rb following adiabatic lattice 
loading with strong harmonic confinement. Solid green line is 
weakly interacting gas and blue triangles are ~3000 lattice- 
confined Rb. (Upper inset) Double occupancy at T/J A = 1 
(filled upside-down triangles) and T/J A = 15 (crosses) for the 
strong harmonic confinement. (Right panels) Integrated den¬ 
sity n(x,y) at points marked I and II at T/J\ = 1 (upper) 
and T/J\ = 15 (lower) for strong harmonic confinement. 


atom number, N A , after the adiabatic loading procedure, 
with filled circles (triangles) for two different trapping 
conditions. The filling at the center of the trap is dis¬ 
played with dashed (solid) lines for the weaker (stronger) 
harmonic confinement at T/J\ = 15. The empty dia¬ 
monds are the experimental results [1]. The rightmost 
panels of Fig. 2 give examples of the resulting distribu¬ 
tions integrated along the z- axis, n(x,y), for N A ~ 4000 
(labeled I) and N A ~ 8000 (labeled II). For these two val¬ 
ues of N A we show the distribution at T/J A = 1 (top), 
which is the low-temperature result, followed by the in¬ 
tegrated density distribution at T/J A = 15 (bottom), 
which is closer to the experimental temperature. 

It is worth noting that while the density at the cen¬ 
ter of the trap is strongly dependent on temperature, f p 
does not display a strong dependence on T/J A within 
the experimentally-relevant temperature range. How¬ 
ever, this does not mean that the efficiency of the for- 
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mation of molecules is unaffected by temperature. This 
is evident if one probes the number of doubly occupied 
sites with increasing temperature. Sites with double (or 
higher) occupancies do not result in molecule formation. 
The upper inset of Fig. 2 shows the double occupancy at 
T/Ja = 1 (inverted triangles) and T/Ja = 15 (crosses). 
In the low-temperature regime ( T/Ja = 1) the first Mott 
shell extends to Na ~ 5000 particles, after which a su¬ 
perfluid region forms at the center of the trap, before 
transitioning to the second Mott shell at Na ~ 8000 (see 
the dashed blue line in Fig 2). However close to or above 
T/Ja = 15, 5 — 10 % of the A atoms are in doubly oc¬ 
cupied sites and hence will not participate in molecule 
formation. 

Next, we consider the fermionic species B of the mix¬ 
ture, taking 40 K for comparison with recent experimental 
data [1]. However, as before, our conclusions are valid for 
generic fermionic species that can be described by a tight- 
binding model. The difference in the polarizability with 
respect to 87 Rb means that the 40 K atoms feel a lattice 
potential of depth Vo = 9 F?b, where Eb is the B species 
recoil energy. At low temperatures, p-wave interactions 
in the spin-polarized gas can be neglected, and the de¬ 
termination of the density reduces to a non-interacting 
problem. As the harmonic trap is separable, this problem 
can be straightforwardly treated by direct diagonaliza- 
tion for the single-particle eigenenergies E n ^ and corre¬ 
sponding wavefunctions | With these single-particle 
quantities we can evaluate the grand canonical partition 
function, from which we find the entropy per particle 
S/N, as well as the on-site density n(r) = n(x)n(y)n(z), 

n (0 = l+exp(/3( 1 £; ri ^-/i)) I^S; (OP 7 ( 3 ) 

where £ = x,y,z and /3 = 1/T is the inverse temperature. 

In Fig. 3 we show the peak filling of species B as a func¬ 
tion of atom number. For fermions we use a Gaussian 
fit with n(x,y ) = \f2Tw~ z f p exp [-x 2 /(2a 2 x ) ~y 2 /(2a 2 )\ 
after integration along z. Here, we also account for the 
experimental imaging resolution and pixelation by apply¬ 
ing a Gaussian filter. The band, delimited by the blue 
triangles and orange squares, shows the range of peak 
fillings for 1 < T^/Jb < 5, with Jb the tunneling of the 
B atoms. As the temperature of the B atoms in the lat¬ 
tice increases, so does the width of the cloud, resulting 
in a decrease in the peak filling. The empty diamonds 
show the experimental data [1]. The plateau indicates 
the atoms forming an incompressible band insulator. On 
the right panel we show the integrated density distribu¬ 
tion at points marked I and II. At II, the B atoms form 
a band insulator in a large region of the lattice. 

Lattice-confined two-species mixture We use a two- 
component mixture of soft-core (“A”) and hard-core 
(“B”) bosons to study the combined effect of interspecies 
interactions and finite temperature on the density of 



FIG. 3. (Color online) Peak filling of the fermionic species B 
vs. atom number. (Main panel) The peak filling of 40 K, used 
as an example for the fermionic species B, at Tq.h./Jb = 1 
(filled squares) and T^/Jb — 5 (filled triangles) extracted 
from the integrated density distributions n(x,y), including 
mimicked resolution and pixelation effects. Results are for 
harmonic trap frequencies (40 x 40 x 260) Hz and a 9Eb deep 
lattice. Empty diamonds show experimental results [1] . Insets 
show the integrated density distributions at points I and II . 

the A species. Hard-core bosons act as a stand-in for 
fermionic B species, as path-integral QMC cannot sim¬ 
ulate fermions due to the sign problem. While there is 
no direct mapping from hard-core bosons to fermions for 
local observables in three dimensions, we compared the 
hard-core and fermionic profiles for the single-species case 
and found excellent agreement. Based on this, we expect 
our results for the local density in the mixture to also be 
valid. 

The two-component mixture is described by the Hamil¬ 
tonian 

H ^ ^ J^{cL^^CLj^ T h.c.) T — 'y ^ ^h,A(^h,A 1) 

(ij), 7 i 

T Fab ^ ^ ^ ^ AG,• (4) 

i i, 7 

Here 7 =A, B and Uab is the interspecies on-site inter¬ 
action which can be tuned to be repulsive or attractive. 
To study the effect of interspecies interactions we have 
performed simulations with —40 < Uab/Ja < 40 with 
10500 A atoms and 25000 B atoms. Previous studies 
have shown that at a given T/ Ja , the presence of attrac¬ 
tive interactions enhances the conversion efficiency [23] . 
This corresponds to an enhancement of f p measured at 
the same T/Ja for every Uab < 0 . In Fig. 4 we show 
f p {T/J A = 0.2) and f p (T/J a = 10) as a function of 
Uab/Ja using filled squares and filled circles connected 
by a line, respectively. We normalize the values to 
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FIG. 4. (Color online) Dependence of species A peak fill¬ 
ing on Uab . (a) Peak filling normalized to zero-temperature 
value, f P (T)/fp , vs. Uab/Ja at T/Ja=0.2 (empty squares) 
and T/Ja—10 (empty circles), ignoring the interaction de¬ 
pendence of the final temperature due to adiabatic loading. 
Accounting for interaction effects during loading at initial 
temperatures E/Tp = 0.1 (0.3) are given with filled squares 
(triangles), (b) The entropy of species A as a function of 
T/Ja is shown for Uab/Ja = 40 and -40 using circles and 
squares, respectively. The dashed line is the entropy of non¬ 
interacting fermions in a harmonic trap, which we use as an 
estimate for Si. (c) Final temperature vs. initial temperature 
for Uab/Ja = —40, -10, 10, and 40 (solid, dash, dot-dash, and 
dot-dot-dash, respectively). We use a harmonic confinement 
of cj r (v z ) = 2tt23(150) Hz 


the peak filling of A at T/Ja = 0.2 in the absence of 
species B. It is clear that attractive interactions lead to 
an enhancement of f p if T/ Ja is kept constant. 

However, this simple analysis does not describe cur¬ 
rent experiments, since starting at the same initial tem¬ 
perature, Ti/Tp , Tp the Fermi temperature, leads to a 
final temperature T^ u /Ja after adiabatic loading that 
depends on the sign and magnitude of Uab- In Fig- 4 
(c) we show Tfi n /JA for Uab/Ja = ~ 40, -10, 10, and 
40 using solid, dash, dot-dash, and dot-dot-dash lines, 
respectively. From Fig. 4(c) it is clear that attractive 
interactions tend to cause more severe heating during 
the adiabatic loading. For the range of initial entropies 
considered in this work, this additional heating inhibits 
and counteracts the benefit gained by making the two 
species attract. In Fig. 4 (a) we show f p as a function of 
Uab/Ja for Ti/Tp =0.1 and 0.3 using filled squares and 
filled triangles, respectively. Following adiabatic loading, 
both attractive and repulsive interactions destabilize the 
A species Mott insulator and reduce the peak filling. 

STIRAP adiabaticity and higher-band transfer. Pairs 
of A and B atoms in an optical lattice site are transferred 
to weakly bound Feshbach molecules following a mag¬ 
neto [11]- or photo [18]-association step. These Feshbach 
molecules are converted into ground state molecules via a 
two-photon STIRAP sequence involving an intermediate 


electronically excited molecular state. As the STIRAP 
process must remove ^100 THz of molecular energy, the 
difference of the wavevectors of the two photons involved 
in the STIRAP sequence, denoted by k u and can have 
significant variation on the lattice scale of a few microns. 
The resulting momentum transfer can excite the resulting 
ground state molecules to higher lattice bands. Here, we 
investigate the adiabaticity of the STIRAP procedure, as 
well as the rate of transfer of molecules to higher bands. 

We consider a STIRAP process which couples the Fesh¬ 
bach molecule (FBM) (|/)) and ground molecular (GSM) 

| g) states through an excited level (|e)) [27]. We assume 
that these are the only states involved and neglect the 
population of other molecular levels or atomic scatter¬ 
ing states during the STIRAP. To match current ex¬ 
periments, we study the case of vanishing single- and 
two-photon detunings. Using adiabatic perturbation the¬ 
ory [28-30] we find that the momentum transfer following 
the process \f) -A \g) has the same form as the case of 
large single-photon detuning, oc e l ^ ku ~ kd )' r . Typical STI¬ 
RAP linewidths are ~ 200 kHz, significantly greater than 
the ~ 20kHz spacing between bands, and so the band 
structure is fully unresolved. Thus, the relative popula¬ 
tion of molecules in the different bands following STIRAP 
are determined by ratios of Rabi frequencies in the basis 
of Wannier states with band index n. We use w n (f) and 
w n {r) to denote the Wannier states for FBM and GSM, 
respectively. For FBM in the lowest band, the relative 
population of GSM in the first excited band along the 
direction £ can thus be estimated by 

(tDi|e tfc ^| Wo ) 2 ^ 2 (Agfa 

{w 0 \e ik ^\w 0 ) ~ n 2 (l+a) 2 ^/V/Pi ’ 

where we have used a harmonic approximation for the 
Wannier functions, with k % = (k u — kg) • £ the momen¬ 
tum transfer along direction £, a the lattice spacing, a 
the polarizability ratio of the GSM to the FBM, V the 
lattice depth for the FBM, and Ep the molecular recoil 
energy. For the parameters of the JILA experiment [1], 
a = 532 nm, k u = 27r/(968 nm) and kd = 27r/(689 nm) 
co-propagating at a 45° angle with respect to the x and 
y lattice axes, and a ~ 0.9, we find the total population 
in the first excited bands to be ^ 1%. In general, our 
results indicate that for experiments in the Lamb-Dicke 
regime k^a/iy/Ep) 1 ^ 1, STIRAP does not induce 
appreciable population transfer to excited bands. 

In conclusion, we studied the combined effects of in¬ 
terspecies interactions, temperature, and adiabatic load¬ 
ing on the successful preparation of a low entropy gas of 
polar molecules. We have shown that interspecies inter¬ 
actions have a significant effect on the final temperature 
of the lattice gas following adiabatic loading, which in 
turn can lead to the depletion of the peak filling of the 
bosonic species and a lower efficiency of molecule for¬ 
mation. Based on our results, molecule formation efH- 
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ciency is greatest when the lattice loading is performed 
with a non-interacting mixture. For lower initial tem¬ 
peratures where the final temperature after loading is 
Tfin/ JA < 5, attractive interactions may lead to a higher 
rate of molecule formation; this temperature is around 
T/Tp < 0.05 for the parameters considered in this work. 
Moreover, we considered the role of the STIRAP pro¬ 
cedure used to convert the resulting Feshbach molecules 
to ground state molecules in the conversion efficiency. 
We find that the probability of promoting molecules to 
higher bands during STIRAP is not a significant concern 
provided the Lamb-Dicke criterion is satisfied. 
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